1 Introduction

The purpose of this document is to provide a reproducible record of all analyses and figures in the main article. The main article is focused on the effect of species’ intrinsic responses on community stability in fluctuating environments. We are going to look at the effect of the distribution of species responses, richness, temperature and nutrients on community temporal stability. Specifically, we are going to look at the effect of fundamental balance (our measurement of the distribution of species’ intrinsic responses) on temporal stability. Species’ intrinsic responses to environmental change can stabilise community biomass in two ways: through response diversity and /or through population stability. Thus, as response diversity is thought to stabilize temporal stability of aggregate community properties via asynchrony, we are going to look at the relationship between response diversity and asynchrony. Subsequently, we are going to look at the relationship between population stability and temporal stability of community biomass. Finally, we use a structural equation model to test the direct and indirect effects of balance on temporal stability of community biomass.

In this document, we also analyse the predictive power of balance and divergence on temporal stability, and we compare the unique explanatory power of balance and divergence. We also look at the interaction between divergence and richness, and between balance and richness. Finally, we assess variable importance using the relative importance of predictors in the full model. This part of the document (section 5), is not included in the main article, but it is an important part of the analysis that justify adopting balance instead of divergence, which was the original metric used to design the experiment.

This document is produced by an Rmarkdown file that includes code to reproduce from data all results presented in the main article.

2 Load datasets, Data wrangling and balance calculation

divergence_df <- read_csv("Data/divergence_df.csv")
load("Data/dens_biomass_poly.RData")

dd_all_pred<-read.csv("Data/morph_dd_pred.csv")
dd_all_pred_nonoise<-read.csv("Data/morph_dd_pred_nonoise.csv")

load("Data/ciliate_traits.Rdata")

df_slopes <- read_csv("Data/df_slopes_cor.csv")

# needs to have id_new variable
ciliate_traits <- ciliate_traits %>%
  dplyr::mutate(
    # Remove dots from the date
    cleaned_date = gsub("\\.", "", date),
    # Extract the part of id after the underscore
    id_suffix = sub(".*_(.*)", "\\1", id),
    # Combine cleaned_date, id_suffix, and species_initial into a new variable
    id_new = paste0(cleaned_date, id_suffix, composition)
  ) %>%
  # Optionally, remove the intermediate columns to clean up
  dplyr::select(-cleaned_date, -id_suffix,-new_id)

uniqueN(ciliate_traits$id_new)==nrow(ciliate_traits) # all unique  ;)

id_dd<-full_join(dd_all_pred,dplyr::select(ciliate_traits,id_new,biomass),join_by("id_new"))


## add day variable

#create a day variable from the date variable

id_dd<-dplyr::mutate(id_dd,date=as.Date(date,format = "%d.%m.%y"))

earliest_date<-min(id_dd$date)
days_since_earliest<-as.numeric(id_dd$date-earliest_date)+1
id_dd<-id_dd%>%dplyr::mutate(day=days_since_earliest)

#create a summarised df on microcosm level with each species seperate
# Make sure, that we have n_frames and not N_frames
names(id_dd)[names(id_dd) == "N_frames"] <- "n_frames"

#extrapolation_factor <- 9.301902  # for 16 x magnification 
extrapolation_factor <- 9.828125  # for 25 x magnification 
video_biomass_species <- c( "C", "P", "S","D","L","T")

biomasses <- id_dd %>%
  dplyr::group_by( day,temperature,nutrients,sample_ID,composition,predict_spec) %>% # group  by xxx
  dplyr::summarize(
    biomass = sum(biomass * n_frames, na.rm = TRUE) / (1 * 125) # if not 3 videos corrections is done below with dens_factor
  ) %>%
  dplyr::mutate(
    biomass = biomass * extrapolation_factor,
    )

biomasses<-biomasses%>%dplyr::mutate(biomass=biomass*1000)

dd_ts_id<-biomasses

#fill up missing dates with biomass<-0

fill_dd<-expand.grid(sample_ID=unique(dd_ts_id$sample_ID),day=unique(dd_ts_id$day),predict_spec=unique(dd_ts_id$predict_spec))
complete_ts<-full_join(fill_dd,dd_ts_id,join_by(sample_ID,day,predict_spec))

complete_ts$biomass[is.na(complete_ts$biomass)]<-0
complete_ts<-complete_ts%>%dplyr::mutate(composition=sub("_.*", "", sample_ID))
complete_ts<-complete_ts %>%
  dplyr::mutate(temperature = sapply(strsplit(as.character(sample_ID), "_"), function(x) paste(x[3], x[4], sep = "-")))
complete_ts<- dplyr::mutate(complete_ts,nutrients = gsub(".*Nut(.*?)_.*", "\\1", sample_ID))

# Now remove wrong combinations of composition and predict_spec / predict_spec

complete_ts<- complete_ts %>%
  rowwise() %>%
  dplyr::filter(predict_spec %in% unlist(strsplit(composition, ""))) %>%
  ungroup()  
complete_ts<-dplyr::mutate(complete_ts,temperature=as.character(temperature),
                    nutrients=as.character(nutrients),
                    richness=nchar(composition))

complete_ts<-complete_ts%>%group_by(sample_ID,composition,day)%>%dplyr::mutate(tot_biomass=sum(biomass))
complete_ts<-complete_ts%>%dplyr::mutate(biom_contribution=biomass/tot_biomass)

df_biomass_mod <- complete_ts

complete_ts<-complete_ts%>%dplyr::mutate(temperature=paste0(temperature," °C"),
                                      nutrients=paste0(nutrients," g/L"))


# introduce slopes of 
names(df_slopes)[names(df_slopes)=="species_initial"]<-"predict_spec"

slope_ts<-full_join(dplyr::select(df_slopes,nutrients,predict_spec,temperature,slope),complete_ts)
slope_ts<-slope_ts%>%dplyr::mutate(w_slope=biom_contribution*slope,
                            sign=sign(slope))

slope_ts<-slope_ts%>%group_by(sample_ID,temperature,nutrients,richness,composition,day,tot_biomass)%>%dplyr::summarize(
  sum_w_slopes=abs(sum(w_slope)),
                   mean_abs_slope=mean(abs(slope)),
  sum_abs_slope=sum(abs(slope)),
  abs_sum_slope=abs(sum(slope)),
  symmetry=abs(sum(sign)))


slope_ts<-slope_ts%>%dplyr::mutate(richness=as.factor(richness))


##create new variable where it checks, where the last observation =0 is; with complete_ts
aggr_ts <- slope_ts %>%
  group_by( sample_ID) %>%
  arrange(day) %>%
  mutate(
    # Create a flag for non-zero tot_biomass
    non_zero_biomass = tot_biomass != 0,
    # Find the last non-zero day
    last_non_zero_day = ifelse(any(non_zero_biomass), max(day[non_zero_biomass], na.rm = TRUE), NA),
    # Find the first zero day after the last non-zero day
    first_zero_day = ifelse(
      !is.na(last_non_zero_day),
      min(day[!non_zero_biomass & day > last_non_zero_day], na.rm = TRUE),
      NA
    ),
    # Flag for days after the first zero day
    is_after_first_zero_day = ifelse(!is.na(first_zero_day), day > first_zero_day, FALSE)
  ) %>%
  ungroup()

aggr_ts<-aggr_ts%>%mutate(rep_var=sub("_[^_]+$", "", sample_ID))

biomass_ts<-aggr_ts%>%group_by(day,temperature,nutrients,richness)%>%summarize(tot_biom=mean(tot_biomass),se_tot_biom=sd(tot_biomass)/sqrt(as.numeric(length(tot_biomass))))

3 Biomass

Let’s have a look at the biomass dynamics in the different environmental treatments.

3.0.1 tot biomass plot

Figure 1 : Community total biomass during the experiment in different environmental treatments. Different color represent richness levels.

4 Main Results

We now look at the main results of the experiment. We are going to look first at the effect of richness, temperature and nutrients on community temporal stability. Then, we are going to look at the relationship between divergence (original response diversity metric) and temporal stability. Finally, we are going to look at the relationship between response diversity and temporal stability.

In the whole analysis, we calculated the temporal stability of total community biomass as the inverse of the coefficient of variation (ICV) (i.e. \(\frac{\sigma}{\mu}\)).

4.0.1 Effect of T, N and R

Figure 2: Effects of richness (a), temperature (b), and nutrients (c) on community total biomass temporal stability.

We can see that richness does not have a clear effect on community temporal stability, while stability was higher at lower temperature, and nutrients increased community temporal stability.

4.0.2 Effect of Divergence

We look at the relationship between divergence (our original response diversity metric) and stability

Figure 3: Relationship between Divergence and temporal stability of total community biomass.

Divergence is positively related to temporal stability, suggesting that response diversity promotes stability. However, the relationship between divergence and stability becomes weaker as richness increases. We think that this is due to divergence considering only the responses of the 2 most “responding” species. Thus, when species richness increases, disregarding the responses of the other species in the community except the 2 responding the most makes the relationship between response diversity and stability weaker.

This is why, after running the experiment, we developed another metric to measure response diversity, which we called balance, and that is presented in the main text of the publication. Balance has several desirable features that makes it a more suitable metric than divergence: Independence of richness, higher predictive power, and accounts for the responses of all species in the community (as opposed to divergence that accounts for only the 2 most “responding” species).

Here, we provide extensive evidence of why balance is a better metric to measure response diversity than divergence, and thus justifying focusing the analysis around balance.

5 Comparing Divergence and Balance

5.1 Predictive power of Divergence and Balance

We first compare how well divergence and balance predict stability (predictive power).

5.1.1 Balance

# 

mod1 <- lm(data=complete_aggr,log10(stability)~log10(balance_f))

# Check model assumptions
#check_model(mod1)

5.1.2 Divergence

mod2 <- lm(data=complete_aggr,log10(stability)~(divergence))

# Check model assumptions
#check_model(mod2)

Table 1: Comparison of model performance of divergence and balance as predictors of stability. Model 1 has balance as predictor and model 2 has divergence as predictor.

model	AIC	AICc	BIC	R2	R2_adjusted	RMSE	Sigma
1	-89.27328	-89.17286	-78.79409	0.1917679	0.1884142	0.1510344	0.1516599
2	-55.71579	-55.61538	-45.23661	0.0720796	0.0682293	0.1618316	0.1625017

A model with Balance as predictor performs better than one with divergence as predictor, and it explains more of the variance in stability than divergence.

Moreover, from Figure 3, it looks like divergence declines in performance as richness increases. Let’s test this analytically. To do than we build a linear model having stability as response variable and either log10(balance) or divergence as predictor for each richness level. We then extract the R squared of the models and their standardised estimates. (standardized estimates were calculated centering divergence and balance using the function scale()).

# getting model estimates for each richness level
lm_divergence_richness_E <- complete_aggr %>%
  nest(data = -richness) %>%
  mutate(
    model = map(data, ~ lm(log10(stability) ~ scale(divergence), data = .x)),
    results = map(model, broom::tidy)
  ) %>%
  unnest(results) %>% dplyr::filter(term=="scale(divergence)") 


# getting model R squared for each richness level

lm_divergence_richness_R <- complete_aggr %>%
  nest(data = -richness) %>%
  mutate(
    model = map(data, ~ lm(log10(stability) ~ scale(divergence), data = .x)),
    results = map(model, broom::glance)
  ) %>%
  unnest(results)

# getting model estimatesf or each richness level
lm_balance_richness_E <- complete_aggr %>%
  nest(data = -richness) %>%
  mutate(
    model = map(data, ~ lm(log10(stability) ~ scale(log10(balance_f)), data = .x)),
    results = map(model, broom::tidy)
  ) %>%
  unnest(results) %>% dplyr::filter(term=="scale(log10(balance_f))") 



# getting model R squared for each richness level
lm_balance_richness_R <- complete_aggr %>%
  nest(data = -richness) %>%
  mutate(
    model = map(data, ~ lm(log10(stability) ~ scale(log10(balance_f)), data = .x)),
    results = map(model, broom::glance)
  ) %>%
  unnest(results)

Figure 4: Performance comparison of divergence vs balance. In (a), the R squared of linear models for divergence and balance are shown for each richness level. In (b), the estimates of the linear models for divergence and balance are shown for each richness level.

We can see that the R squared of divergence as predictor of stability becomes smaller as richness increases, while the R squared of balance as predictor of stability does not (actually increases slightly).

5.2 Comparing unique explanatory power of balance and divergence

Now we build a linear model were stability is modeled as a function of balance and divergence. Then, we compared the variance explained by the full model compared to a model containing either only balance or only divergence.

5.2.1 Full model - balance and divergence

lm_div_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f)+divergence)

# Check model assumptions
# check_model(lm_div_balance)

5.2.2 model with only divergence

lm_div <- lm(data=complete_aggr,log10(stability)~divergence)

# Check model assumptions
# check_model(lm_div)

5.2.3 model with only balance

lm_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f))

# Check model assumptions
# check_model(lm_balance)

5.2.4 Comparision full model vs divergence only and balance only

Table 2: Comparison of model performance of divergence, balance and both as predictors of stability. Model 1 has both balance and divergence as predictors, model 2 has divergence as predictor, and model 3 has balance as predictor.

model	AIC	AICc	BIC	R2	R2_adjusted	RMSE	Sigma
1	-88.97683	-88.80876	-75.00458	0.1974141	0.1907259	0.1505060	0.1514437
1	-55.71579	-55.61538	-45.23661	0.0720796	0.0682293	0.1618316	0.1625017
2	-89.27328	-89.17286	-78.79409	0.1917679	0.1884142	0.1510344	0.1516599

5.2.5 Comparision full model vs balance only

Table 3: Anova table: a model with both balance and divergence as predictors is not significantly different from a model with only balance as predictor.

anova1 <- anova(lm_div_balance,  lm_balance)

# Convert to tidy format
anova_tidy1 <- broom::tidy(anova1)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy1 %>%
  gt() %>%
  cols_label(
    term = "Term",
    sumsq = "Sum of Squares",
    df = "DF",
    statistic = "F Statistic",
    p.value = "p-value"
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      columns = vars(p.value),
      rows = p.value < 0.05
    )
  ) %>%
  tab_options(
    table.width = px(800),            # Adjust table width (e.g., 400px)
    table.font.size = px(12),        # Adjust font size (e.g., 12px)
    data_row.padding = px(10)         # Adjust row padding (e.g., 4px for more compact rows)
  )

Term	df.residual	rss	DF	Sum of Squares	F Statistic	p-value
log10(stability) ~ log10(balance_f) + divergence	240	5.504447	NA	NA	NA	NA
log10(stability) ~ log10(balance_f)	241	5.543171	-1	-0.03872444	1.688429	0.195

5.2.6 Comparision full model vs divergence only and divergence only

Table 4: Anova table: a model with both balance and divergence as predictors is significantly better from a model with only divergence as predictor.

anova2 <- anova(lm_div_balance,  lm_div)


anova_tidy2 <- broom::tidy(anova2)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy2 %>%
  gt() %>%
  cols_label(
    term = "Term",
    sumsq = "Sum of Squares",
    df = "DF",
    statistic = "F Statistic",
    p.value = "p-value"
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      columns = vars(p.value),
      rows = p.value < 0.05
    )
  ) %>%
  tab_options(
    table.width = px(800),            # Adjust table width (e.g., 400px)
    table.font.size = px(12),        # Adjust font size (e.g., 12px)
    data_row.padding = px(10)         # Adjust row padding (e.g., 4px for more compact rows)
  )

Term	df.residual	rss	DF	Sum of Squares	F Statistic	p-value
log10(stability) ~ log10(balance_f) + divergence	240	5.504447	NA	NA	NA	NA
log10(stability) ~ divergence	241	6.364040	-1	-0.8595933	37.47922	0.000

Overall, balance explains more of the variance in stability than divergence, and there is virtually no difference between a model containing only balance and the full model.

5.3 Interaction divergence and richness

Richness had to be transformed to numeric and to be centered to avoid collinearity with divergence

lm_rich_div <- lm(data=complete_aggr,log10(stability)~divergence*scale(as.numeric(richness)))

# check model assumptions
# check_model(lm_rich_div)

Table 5: Type III anova table of the model with divergence and richness as predictors of stability.

anova3 <- car::Anova(lm_rich_div, type = "III")

anova_tidy3 <- broom::tidy(anova3)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy3 %>%
  gt() %>%
  cols_label(
    term = "Term",
    sumsq = "Sum of Squares",
    df = "DF",
    statistic = "F Statistic",
    p.value = "p-value"
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      columns = vars(p.value),
      rows = p.value < 0.05
    )
  ) %>%
  tab_options(
    table.width = px(800),            # Adjust table width (e.g., 400px)
    table.font.size = px(12),        # Adjust font size (e.g., 12px)
    data_row.padding = px(10)         # Adjust row padding (e.g., 4px for more compact rows)
  )

Term	Sum of Squares	DF	F Statistic	p-value
(Intercept)	11.033652088	1	442.55571734	0.000
divergence	0.807044347	1	32.37025122	0.000
scale(as.numeric(richness))	0.001236238	1	0.04958505	0.824
divergence:scale(as.numeric(richness))	0.249582101	1	10.01064605	0.002
Residuals	5.958668583	239	NA	NA

Divergence significantly interact with richness, suggesting that the relationship between divergence and stability changes with richness. While an ideal metric of response diversity should be independent of richness.

We repeat the same model using balance instead of divergence.

lm_rich_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f)*scale(as.numeric(richness)))

# check model assumptions
# check_model(lm_rich_balance)

Table 6: Type III anova table of the model with balance and richness as predictors of stability.

anova4 <- car::Anova(lm_rich_balance, type = "III")

anova_tidy4 <- broom::tidy(anova4)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy4 %>%
  gt() %>%
  cols_label(
    term = "Term",
    sumsq = "Sum of Squares",
    df = "DF",
    statistic = "F Statistic",
    p.value = "p-value"
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      columns = vars(p.value),
      rows = p.value < 0.05
    )
  ) %>%
  tab_options(
    table.width = px(800),            # Adjust table width (e.g., 400px)
    table.font.size = px(12),        # Adjust font size (e.g., 12px)
    data_row.padding = px(10)         # Adjust row padding (e.g., 4px for more compact rows)
  )

Term	Sum of Squares	DF	F Statistic	p-value
(Intercept)	9.54928736	1	414.779284	0.000
log10(balance_f)	1.26534818	1	54.961191	0.000
scale(as.numeric(richness))	0.02471247	1	1.073401	0.301
log10(balance_f):scale(as.numeric(richness))	0.04049552	1	1.758948	0.186
Residuals	5.50239554	239	NA	NA

Balance does not significantly interact with richness, suggesting that the relationship between balance and stability is stable across richness levels.

5.4 Variable importance

Finally, we assess variable importance using the relative importance of predictors in the full model. We use the package vip (https://cran.r-project.org/web/packages/vip/vignettes/vip.html) to calculate the relative importance of predictors in the full model. The function vip::vip for multiple linear regression, or linear models (LMs), uses the absolute value of the -statistic as a measure of VI. Motivation for the use of the associated 𝑡-statistic is given in Bring (1994) [https://www.tandfonline.com/doi/abs/10.1080/00031305.1994.10476059].

vip::vip(lm_div_balance)

Figure 5: Variable importance in the model including both balance and divergence as predictors of stability.

We believe that the extensive evidence here provided justifies focusing the analysis around balance, and not divergence, as a metric of response diversity. We will thus only look at balance for the rest of the analysis.

6 Effect RD

We are now going to look at how response diversity (balance) affected temporal stability of total community biomass. We are going to look at the relationship between fundamental balance (so based only on species response surfaces measured in monoculture), an realised balance (measured accounting for species contribution to balance).

This is fundamentally testing our most important hypothesis.

Figure 6: Effects of fundamental and realised response diversity (measured as balance) on total community biomass temporal stability.

We can see that balance is always negatively related to temporal stability, which means that response diversity promotes stability across richness levels. Interestingly, we see that there is little difference between fundamental and realised balance. Yet, as the richness increases, the relationship between realised balance and stability becomes steeper compared to fundamental balance.

But is the difference in the slope of fundamental and realised balance significant? We can test this using a linear model with interaction between balance and type (factor: fundamental or realised balance).

6.1 Balance: realised vs fundamental

# compare if the slope of fundamental and realised balance is significantly different for each richness level
# Fit the linear model with interaction
model_F <- lm(log10(1/CV) ~ log10(balance_f) , data = complete_aggr)
model_R_F <- lm(log10(1/CV) ~ log10(balance_f) +log10(balance_r) , data = complete_aggr)

Table 7: Anova table of the model with only realised balance vs one with both realised and fundamental balance as predictors of stability.

anova5 <- anova(model_F, model_R_F)

anova_tidy5 <- broom::tidy(anova5)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy5 %>%
  gt() %>%
  cols_label(
    term = "Term",
    sumsq = "Sum of Squares",
    df = "DF",
    statistic = "F Statistic",
    p.value = "p-value"
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      columns = vars(p.value),
      rows = p.value < 0.05
    )
  ) %>%
  tab_options(
    table.width = px(800),            # Adjust table width (e.g., 400px)
    table.font.size = px(12),        # Adjust font size (e.g., 12px)
    data_row.padding = px(10)         # Adjust row padding (e.g., 4px for more compact rows)
  )

Term	df.residual	rss	DF	Sum of Squares	F Statistic	p-value
log10(1/CV) ~ log10(balance_f)	241	5.543171	NA	NA	NA	NA
log10(1/CV) ~ log10(balance_f) + log10(balance_r)	240	5.522059	1	0.02111215	0.9175772	0.339

Table 7: Linear model results for the interaction between balance and type (fundamental vs realised) as predictors of stability for richness level.

model	AIC	AICc	BIC	R2	R2_adjusted	RMSE	Sigma
1	453.8407	453.9411	464.3198	0.1917679	0.1884142	0.1510344	0.1516599
2	454.9134	455.0814	468.8856	0.1948462	0.1881365	0.1507466	0.1516858

A model with both fundamental and realised balance as predictors improved very little the variance explained by the model.

7 Linear models

7.1 Model: Fundamental balance and Interaction between temperature and nutrients

We may expect and interactive effect of the environmental variables on stability. We thus build a linear model with interaction between temperature and nutrients. However, there is high collinearity between temperature and nutrients, which may affect the model results.

So we transformed nutrients and temperature to numeric, and transformed temperature regimes in values = 1, 2, 3. Then, we centered the variables to avoid collinearity with the interaction term.

# transform nutrients and temperature to numeric. For this the units need to be removed, and temperature regimes should be transformed in values = 1, 2, 3
complete_aggr_2<- complete_aggr %>%
  # Remove the units from the 'nutrients' and 'temperature' columns
  mutate(
    nutrients = as.numeric(gsub(" g/L", "", nutrients)),  # Convert nutrients to numeric
    temperature = gsub(" °C", "", temperature)            # Remove the unit but keep as character
  ) %>%
  # Convert temperature ranges to numeric codes using case_when
  mutate(
    temperature = case_when(
      temperature == "18-21" ~ 1,
      temperature == "22-25" ~ 2,
      temperature == "25-28" ~ 3,
      TRUE ~ NA_real_         # Handle unexpected values with NA
    )
  )


# Fit the linear model with interaction
lm_full_int<-lm(data=complete_aggr_2,log10(stability)~log10(balance_f)+scale(nutrients)*scale(temperature)+richness)

# check model assumptions
# check_model(lm_full_int)

#summary(lm_full_int)

Table 8: Linear model results for the effects of balance, richness, nutrients, and temperature on community stability. Estimates are presented with 95% confidence intervals and p-values.

Predictor	Linear Regression Results	95% CI¹	Linear Regression Results
Predictor	Estimate	95% CI¹	p-value
log10(balance_f)	-0.03	-0.07, 0.00	0.043
scale(nutrients)	0.09	0.07, 0.10	<0.001
scale(temperature)	-0.05	-0.07, -0.03	<0.001
richness
richness3 - richness2	-0.04	-0.09, 0.00	0.082
richness4 - richness2	-0.02	-0.06, 0.03	0.7
richness4 - richness3	0.03	-0.02, 0.07	0.4
scale(nutrients) * scale(temperature)	-0.01	-0.03, 0.01	0.2
¹ CI = Confidence Interval

The relationship between community stability and the predictors, including balance, nutrient and temperature levels, and species richness, was analyzed using a linear model (Table 1). Overall, the model explained 47.2% of the variation in community stability (adjusted R2 = 0.4581, F6,236 = 35.09, p<2.2×10−16). The results showed that the balance of species’ responses to environmental conditions (log10(balancef)had a small but significant negative effect on stability.

Nutrient availability had a strong, positive effect on community stability (β=0.088±0.0088, p<2×10−16p. In contrast, temperature significantly reduced stability (β=−0.049±0.0100.010β=−0.049±0.010, p=6.14×10−6). Species richness had mixed effects on stability. Communities with richness of three species exhibited significantly lower stability compared to the reference group (β=−0.042±0.019, p=0.0323). However, richness of four species did not significantly affect stability (β=−0.016±0.020, p=0.416). The interaction between nutrients and temperature was not statistically significant (β=−0.011±0.009, p=0.195), suggesting that their combined effects on stability were negligible under the tested conditions.

Table 9: ANOVA table of the model with interaction between temperature and nutrients as predictors of stability.

Term	DF	Sum of Squares	meansq	F Statistic	p-value
ANOVA Table for Linear Model
log10(balance_f)	1	1.31521873	1.31521873	85.633282	0.000
scale(nutrients)	1	1.50988193	1.50988193	98.307712	0.000
scale(temperature)	1	0.30972717	0.30972717	20.166192	0.000
richness	2	0.07294304	0.03647152	2.374644	0.095
scale(nutrients):scale(temperature)	1	0.02595802	0.02595802	1.690115	0.195
Residuals	236	3.62466104	0.01535873	NA	NA

8 Asynchrony

Response diversity (aka balance) has been suggested as a mechanism that promotes temporal stability of community biomass by promoting species asynchrony.

We thus calculated the asynchrony index suggested by Gross et al. 2014 to calculate the effect of asynchrony on temporal stability and to see how reponse diversity relate to asynchrony. The index ranges between -1 and 1, with -1 indicating perfect asyncrony and 1 being perfectly synchronous, and 0 indicating random variation.

8.0.1 Plot stability vs. Asynchrony Gross

Figure 8: Relationship between temporal stability and asynchrony (Gross) divided by nutrient level.

The Pearson’s correlation between asynchrony and stability is significant (estimate = -0.23, p < 0.001).

cor.test((-1*async_aggr$synchrony_Gross),async_aggr$stability)

## 
##  Pearson's product-moment correlation
## 
## data:  (-1 * async_aggr$synchrony_Gross) and async_aggr$stability
## t = 3.7927, df = 239, p-value = 0.0001888
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1153693 0.3539711
## sample estimates:
##       cor 
## 0.2382622

8.0.2 Plot Asynchrony Gross vs fundamental balance

Figure 9: Relationship between asynchrony (Gross) and fundamental balance divided by nutrient level.

The Pearson’s correlation between asynchrony and balance is significant (estimate = 18, p = 0.003).

cor.test((-1*async_aggr$synchrony_Gross),(async_aggr$balance_f))

## 
##  Pearson's product-moment correlation
## 
## data:  (-1 * async_aggr$synchrony_Gross) and (async_aggr$balance_f)
## t = -2.9796, df = 239, p-value = 0.003184
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3082462 -0.0644258
## sample estimates:
##        cor 
## -0.1892515

9 Population stability

The relationship between community stability and the stability of the individual populations that make up the community is a key question in community ecology. Importantly, community stability can result from low population stability, if populations fluctuate asynchronously, or from high population stability, if populations do not fluctuate much. Synthesis of the literature suggests diversity can have a positive or negantive effect on population stability Campbell et al 2010 and (Xu et al 2021)[https://onlinelibrary.wiley.com/doi/full/10.1111/ele.13777].

Theoretical work has suggested that community stability is a product of two quantities: the (a)synchrony of population fluctuations, and an average species-level population stability that is weighted by relative abundance Thibaut & Connolly 2013.

Critically, a balance value close to zero can result from high response diversity, but also from high population stability (population biomass does not change largely over time). We want to look now at whether our new metric of balance can capture these two stabilising mechanisms.

Thus, we first calculate species-level population stability weighted by relative abundance.

–>

10 SEM

Finally, we use a structural equation model (SEM) to explore how stability is influenced by asynchrony, population stability, balance and, nutrient levels. In order to develop a hypothesis regarding the influence of stability, we have drawn on existing literature. This has enabled us to posit that stability is influenced by two key factors: asynchrony and population stability. In turn, these are influenced by balance and, in our particular case, by nutrient levels.

## lavaan 0.6-19 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        12
## 
##   Number of observations                           241
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                                 2.502       2.362
##   Degrees of freedom                                 3           3
##   P-value (Chi-square)                           0.475       0.501
##   Scaling correction factor                                  1.059
##     Satorra-Bentler correction                                    
## 
## Model Test Baseline Model:
## 
##   Test statistic                               852.984     819.012
##   Degrees of freedom                                 9           9
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  1.041
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    1.000       1.000
##   Tucker-Lewis Index (TLI)                       1.002       1.002
##                                                                   
##   Robust Comparative Fit Index (CFI)                         1.000
##   Robust Tucker-Lewis Index (TLI)                            1.002
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)                516.122     516.122
##   Loglikelihood unrestricted model (H1)             NA          NA
##                                                                   
##   Akaike (AIC)                               -1008.245   -1008.245
##   Bayesian (BIC)                              -966.427    -966.427
##   Sample-size adjusted Bayesian (SABIC)      -1004.464   -1004.464
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.000       0.000
##   90 Percent confidence interval - lower         0.000       0.000
##   90 Percent confidence interval - upper         0.101       0.097
##   P-value H_0: RMSEA <= 0.050                    0.692       0.723
##   P-value H_0: RMSEA >= 0.080                    0.127       0.107
##                                                                   
##   Robust RMSEA                                               0.000
##   90 Percent confidence interval - lower                     0.000
##   90 Percent confidence interval - upper                     0.102
##   P-value H_0: Robust RMSEA <= 0.050                         0.703
##   P-value H_0: Robust RMSEA >= 0.080                         0.126
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.018       0.018
## 
## Parameter Estimates:
## 
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                      Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   stability ~                                                             
##     asynchrny_Grss      0.157    0.011   14.493    0.000    0.157    0.326
##     pop_stability       0.957    0.022   42.596    0.000    0.957    0.987
##   asynchrony_Gross ~                                                      
##     log_balance_f      -0.093    0.034   -2.767    0.006   -0.093   -0.185
##     nutrients          -0.211    0.023   -9.072    0.000   -0.211   -0.503
##   pop_stability ~                                                         
##     log_balance_f      -0.073    0.011   -6.852    0.000   -0.073   -0.292
##     nutrients           0.136    0.008   17.854    0.000    0.136    0.653
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .stability         0.098    0.011    8.847    0.000    0.098    0.594
##    .asynchrny_Grss   -0.090    0.055   -1.625    0.104   -0.090   -0.262
##    .pop_stability    -0.630    0.018  -34.512    0.000   -0.630   -3.702
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .stability         0.003    0.000    7.632    0.000    0.003    0.099
##    .asynchrny_Grss    0.088    0.011    8.214    0.000    0.088    0.753
##    .pop_stability     0.012    0.001   11.568    0.000    0.012    0.406
## 
## R-Square:
##                    Estimate
##     stability         0.901
##     asynchrny_Grss    0.247
##     pop_stability     0.594

Figure 10.1: SEM.

Model Fit Indices The model fit indices suggest that the model fits the data well.

Chi-Square Test (User Model): The chi-square test statistic for the user model is χ 2 =1.626 (scaled = 1.465) with 3 degrees of freedom and a p-value of 0.653 (scaled = 0.690). This indicates a good fit, as the test is non-significant, suggesting no significant difference between the observed and model-implied covariance matrices.

Comparative Fit Index (CFI) and Tucker-Lewis Index (TLI): Both CFI and TLI values are 1.000, indicating an excellent model fit. Values close to or above 0.95 are generally considered good.

Root Mean Square Error of Approximation (RMSEA): The RMSEA is 0.000, with a 90% confidence interval ranging from 0 to 0.090 (scaled = 0.080). This indicates a very good fit, as RMSEA values below 0.05 are ideal, and values below 0.08 are acceptable. The p-values for the RMSEA hypothesis tests suggest strong support for a close fit (RMSEA <= 0.05) and little evidence for a poor fit (RMSEA >= 0.08).

Standardized Root Mean Square Residual (SRMR): The SRMR value is 0.017, which is also within the acceptable range (values below 0.08 are generally considered good). Overall, the fit indices suggest that the model is an excellent fit for the data.

Regression Paths and Interpretation

Stability Regressions

Stability ~ Asynchrony_Gross (asynchrny_Grss): The standardized estimate for the effect of asynchrony on stability is 0.340 (p < 0.001), indicating a significant positive association. Higher asynchrony in species dynamics is associated with increased community stability.

Stability ~ Population Stability (pop_stability): The standardized estimate is 0.977 (p < 0.001), showing a strong positive relationship. This suggests that community stability is highly dependent on the stability of individual populations within the community.

Asynchrony_Gross Regressions

Asynchrony_Gross ~ Log10(Balance): The standardized estimate is -0.176 (p = 0.013), indicating a significant negative effect. Higher balance leads to lower asynchrony, suggesting that as balance increases, species within the community fluctuate more synchronously.

Asynchrony_Gross ~ Nutrients: The standardized estimate is -0.469 (p < 0.001), showing a strong negative relationship. Higher nutrient levels appear to reduce asynchrony, possibly by causing similar responses across species.

Population Stability Regressions

Population Stability ~ Log10(Balance): The standardized estimate is -0.296 (p < 0.001), indicating that higher balance is associated with lower population stability.

Population Stability ~ Nutrients: The standardized estimate is 0.635 (p < 0.001), showing that higher nutrient levels are associated with increased population stability, likely because nutrients enhance conditions that support stable population dynamics.

Variances and R-Squared Values R-Squared for Stability: The model explains 90.4% of the variance in community stability, indicating strong predictive power.

R-Squared for Asynchrony_Gross: The model explains 21.9% of the variance in asynchrony, which is moderate.

R-Squared for Population Stability: The model explains 56.2% of the variance in population stability, showing that nutrients and balance are important but not the only factors influencing it.

Summary Interpretation Model Fit: The model has an excellent fit, as indicated by the fit indices. Stability: Community stability is strongly influenced by both population stability and asynchrony among species, with population stability being the stronger predictor. Asynchrony and Balance: Asynchrony decreases with increasing balance and nutrients, suggesting that these factors promote more synchronized fluctuations among species. Population Stability and Nutrients: Higher nutrient levels are associated with increased population stability, suggesting that nutrient availability supports stable population dynamics. Conversely, higher balance is associated with decreased population stability.

Extended results for: The balance of nature: Critical Role of Species Intrinsic Responses for Stability

Til Hämmig, Francesco Polazzo, Owen L. Petchey, Frank Pennekamp

21 November, 2024